function creepplot(filenumber)
count=filenumber;

N=readdata(count,1);
S=readdata(count,2);
T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
time=readdata(count,6);
%     dtnom=readdata(count,7);
dx=readdata(count,8);
dz=readdata(count,9);
mx=readdata(count,10);
mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
mtype=readdata(count,14);
mTemp=readdata(count,15);
mstressxx=readdata(count,16);
mstressxz=readdata(count,17);
mstrainxx=readdata(count,18);
mstrainxz=readdata(count,19);
%mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
P=readdata(count,23);
Nsurfaces=readdata(count,24);
surfaces=readdata(count,25);

time/(60*60*365*24*1e6)
x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end

zP=zeros(1,T)+0.5*dz(1); xP=zeros(1,S)+0.5*dx(1);
for t=1:T-2 zP(t+1)=zP(t)+0.5*(dz(t)+dz(t+1)); end
zP(T)=zP(T-1)+dz(T-1);
for s=1:S-2 xP(s+1)=xP(s)+0.5*(dx(s)+dx(s+1)); end
xP(S)=xP(S-1)+dx(S-1);
[ZP,XP] = meshgrid(zP,xP); P=P(1:end-1,1:end-1); XP=XP(1:end-1,1:end-1); ZP=ZP(1:end-1,1:end-1);
mP=interp2(XP',ZP',P',mx,mz);

figure(1); clf;

mcreepstrain=nan(1,N); mdislstrain=nan(1,N); mdiffstrain=nan(1,N);
for e=1:N
    if (mtype(e)==3 | mtype(e)==4) && ~isnan(mP(e))
        [mcreepstrain(e) mdislstrain(e) mdiffstrain(e)]=serialvisc(3.5,0.5,2.4169e-016,540e3,mP(e),15e-6,mTemp(e),3.3482e+008*exp((300e3+5e-6*mP(e))/(8.314472*(mTemp(e)+273))),mstrainxx(e),mstrainxz(e),mstressxx(e),mstressxz(e));
    end
end
    
subplot(3,1,1);
[grid,xgrid,zgrid]=bilingrid(log10(mcreepstrain),mx,mz,dx,dz,1,1);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 2000 -40 670]); title(['log creep strain rate']); ylabel('Depth (km)');

subplot(3,1,2);
[grid,xgrid,zgrid]=bilingrid(log10(mdislstrain),mx,mz,dx,dz,1,1);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 2000 -40 670]); title(['log dislocation creep strain rate']); ylabel('Depth (km)');

subplot(3,1,3);
[grid,xgrid,zgrid]=bilingrid(log10(mdiffstrain),mx,mz,dx,dz,1,1);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 2000 -40 670]); title(['log diffusion creep strain rate']); ylabel('Depth (km)');

figure(2); clf;
[grid,xgrid,zgrid]=bilingrid(log10(mdislstrain./mdiffstrain),mx,mz,dx,dz,1,1);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 2000 -40 670]); title(['log (dislocation creep rate / diffusion creep rate)']); ylabel('Depth (km)');


function [creepstrain dislstrain diffstrain]=serialvisc(n,F,Ad,E,P,V,T,linvisc,strainxx,strainxz,oldstressxx,oldstressxz)
R=8.314472;
T=T+273;
my=6.7e10;

dt=0.5e4*60*60*24*365;

sigmaII=secinv(oldstressxx,oldstressxz);
lowsigma=sigmaII; highsigma=sigmaII;
visc=(2/sigmaII)^(n-1)*F^n/Ad*exp((E+P*V)/(R*T));
if (linvisc>0) visc=1/(1/visc+1/linvisc); end
[stressxx,stressxz]=calcstress(oldstressxx,oldstressxz,strainxx,strainxz,visc,my,dt);
sigmaII=secinv(stressxx,stressxz); %plot(visc,sigmaII,'r.'); set(gca,'yscale','log','xscale','log');
if sigmaII>highsigma highsigma=sigmaII; end
if sigmaII<lowsigma lowsigma=sigmaII; end
visc=(2/sigmaII)^(n-1)*F^n/Ad*exp((E+P*V)/(R*T));
if (linvisc>0) visc=1/(1/visc+1/linvisc); end
[stressxx,stressxz]=calcstress(oldstressxx,oldstressxz,strainxx,strainxz,visc,my,dt);
sigmaII=secinv(stressxx,stressxz); %plot(visc,sigmaII,'r.');
if sigmaII>highsigma highsigma=sigmaII; end
if sigmaII<lowsigma lowsigma=sigmaII; end



i=0;
while highsigma/lowsigma>1.01;
    i=i+1;
    %sigmaII=(highsigma*lowsigma)^0.5; sigmaold=sigmaII;
    sigmaII=(highsigma+lowsigma)*0.5; sigmaold=sigmaII;
    visc=(2/sigmaII)^(n-1)*F^n/Ad*exp((E+P*V)/(R*T));
    if (linvisc>0) visc=1/(1/visc+1/linvisc); end
    [stressxx,stressxz]=calcstress(oldstressxx,oldstressxz,strainxx,strainxz,visc,my,dt);
    sigmaII=secinv(stressxx,stressxz);
    if sigmaII>sigmaold lowsigma=sigmaold; else highsigma=sigmaold; end
end

creepstrain=sigmaII/visc/2-(sigmaII-secinv(oldstressxx,oldstressxz))/dt/my/2;
diffstrain=sigmaII/linvisc/2;
dislstrain=creepstrain-diffstrain;
end

function [stressxx,stressxz]=calcstress(oldstressxx,oldstressxz,strainxx,strainxz,visc,my,dt)
stressxx=2*visc*strainxx*dt*my/(dt*my+visc)+oldstressxx*(visc/(dt*my+visc));
stressxz=2*visc*strainxz*dt*my/(dt*my+visc)+oldstressxz*visc/(dt*my+visc);
end

function secinv=secinv(stressxx,stressxz)
secinv=(stressxx^2+stressxz^2)^0.5;
end

end

